Linear response strength functions with iterative Arnoldi diagonalization 
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We report on an implementation of a new method to calculate RPA strength functions with 
iterative non-hermitian Arnoldi diagonalization method, which does not explicitly calculate and 
store the RPA matrix. We discuss the treatment of spurious modes, numerical stability, and how the 
method scales as the used model space is enlarged. We perform the particle-hole RPA benchmark 
calculations for double magic nucleus 132 Sn and compare the resulting electromagnetic strength 
functions against those obtained within the standard RPA. 
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I. INTRODUCTION 

The linear response theory (LRT) obtained from the 
linearized time-dependent mean field method is an im- 
portant tool for calculating properties of excited states 
of many-fermion systems, such as nuclear giant reso- 
nances. In its charge-changing version, it can also give 
access to the beta-decay strengths. This method is espe- 
cially important in heavy nuclei, where the shell-model or 
configuration-interaction approaches are intractable. An 
advantage is also that LRT does not require the knowl- 
edge of an interaction and can therefore be used both 
within density functional theory (DFT) and phenomeno- 
logical energy density functional (EDF) approaches, giv- 
ing rise to a set of equations of RPA type. Below, for 
simplicity we refer to this method and associated equa- 
tions simply as RPA method/equations. Strength func- 
tions obtained in this way probe new aspects of the EDFs 
and thus have a potential of constraining parameters in 
phcnomenological nuclear EDFs. 

The purpose of the present study is to present an im- 
plementation of an efficient RPA algorithm that is based 
on the local nuclear EDF. For electronic systems, simi- 
lar methods have been used since many years, see, e.g., 
the recent Ref. [l| for a review, and they also consti- 
tute parts of standardized computer packages such as 
GAMESS @, [|. There are two essential elements of 
these methods, which are at the heart of their efficiency 
and scalability, namely, (i) the RPA equations are solved 
iteratively and (ii) the RPA matrix does not have to be 
explicitly calculated. The second of these elements is par- 
ticularly important; it is based on the observation that 
the action of the RPA matrix on the vector of RPA am- 
plitudes can proceed through the calculation of the mean 
fields corresponding to these amplitudes. 

In nuclear physics context, probably the first study 
that used the concept of mean fields in the RPA method 
was that by P.-G. Reinhard Q. Iterative solutions of 
the RPA equations were introduced by Johnson, Bertsch, 
and Hazelton [f|, and applied to the case of separable 



interactions, but in fact these methods can also be ap- 
plied in more complicated situations, as we show here. 
Strangely enough, these very efficient methods have not 
yet been used in practical applications. Only very re- 
cently, Nakatsukasa et al. [y, 0] have implemented the 
analogous approach within the so-called finite amplitude 
method (FAM). 

Our present implementation pertains to the spherical 
symmetry with neglected pairing correlations - thus it 
constitutes only a proof-of-principle study. The real chal- 
lenge is in solving the quasiparticle RPA (QRPA) prob- 
lem in deformed nuclei. Although at present, a few imple- 
mentations that are based on solving the standar d Q RPA 
equations already exist [8, 9] or begin to emerge [HI HH , 
such a route is bound to be blocked by the shear di- 
mensionality of the problem. On the other hand, as we 
show here, methods based on the iterative solutions using 
mean fields have much better scalability properties and 
are potentially very promising. 

The paper is organized as follows. In Sees. IH1 and IIIII 
we lay down the essential features of the method by pre- 
senting the use of mean fields and iterative solution of the 
RPA method, respectively. Then, in Sec. IIVI we present 
the method to remove the spurious RPA states, which is 
tailored to be used within the iterative approach. Sec- 
tions [V] and I VII present the convergence and scalability 
properties of our method, respectively, and summary and 
conclusions are given in Sec. IVII1 



II. RPA FROM LINEARIZED TDHF 

To be concise, below we present a less general deriva- 
tion than the standard method to derive the RPA 
equations from linearized time-dependent Hartree-Fock 
(TDHF) equations. For density independent forces or 
functionals with terms quadratic in density, the density 
matrix and mean field of a time-dependent nuclear state 
are expressed as 



Pit) 



(1) 



2 



h(t) = ho + h^ + hl 
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where po and are the Hartree-Fock (HF) ground-state 
density matrix and mean field, respectively. Inserting 
p(t) and h{t) of Eqs. © and © into the TDHF equa- 
tion, and keeping only terms linear in the fluctuating 
quantities p and h, we get a linearized TDHF equation, 
or the RPA equations: 
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where we use the letter m for particle states and i for 
hole states, and where e m ^ are the HF single-particle 
energies. The fields h are the first functional derivatives 
of the used EDF, evaluated using the density amplitudes 
p of Eq. |T]). Density dependence of the used EDF beyond 
quadratic gives rise to rearrangement fields in h. These 
rearrangement parts of h must be linearized around po 
to make our RPA equations explicitly first order in p. 

One way to achieve this is by calculating functional 
derivatives of the rearrangement parts of h with respect 
to density, which technically makes our mean-field rou- 
tine differ from the standard HF routines. Since in our 
implementation we use the standard Skyrme forces that 
have simple density dependencies of the coupling con- 
stants, the explicit functional differentiation does not 
cause any mathematical or performance problems. Had 
a EDF with more complex density dependence been used 
it would have been an advantage to instead use the FAM 
method 6] for linearization. 

If the matrix elements of h in Eqs. (O and (HJ) are 
expanded in terms of the particle-hole (p-h) and hole- 
particle (h-p) matrix elements of p, wc obtain the tra- 
ditional RPA equations. In this work, we do not con- 
struct the RPA matrix, but directly solve Eqs. ((3]) and 
(@| by calculating the matrix elements of fields h us- 
ing a HF mean-field routine that uses the time-reversal- 
invariance breaking density matrix p. Since the same 
routine is used to evaluate the HF and RPA mean fields, 
the method is always fully self-consistent [3, [HI]. In 
the following equations, we use the standard abbrevia- 
tions X ml = p Ul mi and Y ml = p u ,im- The density vector 
that contains the p-h matrix elements of p u is defined 
as X u = {X muil ,X m2t i 2 ,. . .,X mD>iD ), and similarly for 
the vector y of h-p elements, where D is the number of 
allowed p-h configurations. Overlaps of RPA vectors are 
defined as 
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(X,Y\X',Y') = ( x*,y* ) 



(5) 



and the minus sign results from the RPA norm matrix. 



III. ITERATIVE SOLUTION OF THE RPA 
EQUATIONS 

The RPA equations ^ and Q constitute a non- 
hermitian eigenproblem with non-positive-definite norm. 



We solve this problem by using an iterative method that 
during each iteration only needs to know the product of 
the RPA matrix and a density vector, that is, the right- 
hand sides of Eqs. (SJ and Q: 



= (e m -e i )X^ i + h mi (^ k ,y k ), 
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where index k labels iterations and the mean fields 
h(X k , y k ) depend linearly on the density vectors X k and 
y k . Expressed through the standard RPA matrices A 
and B pJI], Eqs. ([6]) and ([7]) for a positive norm basis 
vector and for its opposite norm partner vector read: 
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In exact arithmetic A — A' and B = B' and therefore 
either Eqs. ([5]) or ^ could be used in the iteration pro- 
cedure with equivalent results. Nevertheless, below we 
use them both to stabilize the iteration process. 

Various iterative methods, which only need to know 
the products of the diagonalized matrix and vectors, ex- 
ist for non-hermitian matrix eigenvalue equations, and 
good examples with pseudocode are shown in Ref . [l| . We 
chose the non-hermitian Lanczos method [5j in a modified 
form, because it conserves all odd-power energy weighed 
sum rules (EWSR) if the starting vector (pivot) of it- 
eration is chosen correctly. In this work, we start from 
a pivot vector that has its elements set to the matrix 
elements of electromagnetic multipole operator, 



X,, 



t \r p Y JM \cf>i), YL = 0, 



(10) 



where p = 2 and J = for the + mode, p = 1 and J = 1 
for the 1~ mode, and p = 2 and J = 2 for the 2 + mode. 
The constant iV 1 is used to normalize the pivot vector 
to unity. For the IS 1 — mode we only present results 
obtained with the operator: 
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to stay consistent with Refs. [l6[ and [T^- For the choice 
of pivot in (|10p one can prove 0] that all odd-power 
EWSRs are satisfied throughout the iteration procedure. 

Because we calculate the RPA matrix- vector products 
by using the mean-field method, and not with a precal- 
culated RPA matrix, we introduce small, but significant 
numerical noise to the resulting vectors. If corrective 
measures are not used to remove or reduce this noise, 
the iteration method fails and produces complex RPA 
eigenvalues early on in the iteration. We stabilize our 
iterative solution method by modifying the method of 
Ref. Q in two ways. First, we use the non-hermitian 
Arnoldi method instead of the non-hermitian Lanczos 
method. The advantage of Arnoldi method is that it 
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orthogonalizes each new basis vector against all previous 
basis vectors and their opposite norm partners, that is, 
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where the overlap matrices , bik , a'* k , and b'* k are calcu- 
lated as in Eq. ([5]). Again, in exact arithmetic, Eqs. (TT2"]) 
and (|13[) are equivalent and the lower matrix elements a!* k 
and b'* k in Eq. (fT3"|) are exact complex conjugates of the 
elements of the upper Krylov-space [18j RPA matrices. 
We assume this to keep the Krylov-space RPA matrix in 
the standard form. 

In the Lanczos method, only the tridiagonal parts of 
RPA matrices are calculated, and small changes in ba- 
sis vectors due to Lanczos re-orthogonalization (which 
always must be used to preserve orthogonality) do not 
show up in the constructed RPA matrix. In the Arnoldi 
method, these small but important elements outside the 
tridigonal part improve the stability as compared to 
Lanczos. 

The norm of the obtained new residual vector in 
Eq. (fT2"|) can be either positive or negative. We do not 
in practice use Eq. (|13p . which in exact arithmetic would 
duplicate the results of Eq. (fT2"]) . Instead, we store only 
the positive-norm basis states and use a similar method 
as in Ref. Q to change sign of the norm in case the 
norm of the residual vector in Eq. (|12p is negative. Thus, 
explicitly, for the positive norm of the residual vector 
N k+1 = ' (X k+1 ,Y k+1 \X k+1 ,Y k+1 ), we define the new 
normalized positive-norm basis vector as 
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If N k+1 < 0, the new normalized positive norm basis 
vector is defined as 
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When maximum number of iterations has been made 
or the iteration has been stopped, the generated Krylov- 
space RPA matrix, with dimension d « D, is diagonal- 
ized with standard methods, that is we solve: 
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The approximate RPA solutions are then obtained by 
transforming the Krylov-space basis vectors, 
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for all k = 1 , . . . , d, and these vectors are used to evaluate 
the strength functions. 

Standard RPA method that constructs and diagonal- 
izes the full RPA matrix can ensure that the lower matri- 
ces in the RPA supermatrix are exact complex conjugates 
of the upper matrices. Our mean-field method can have 
small differences in the implicitly used upper and lower 
RPA matrices due to finite numerical precision. The con- 
sequence of this is that we will in general have a'^ =/= aij 
and bjA ^ by. This spoils the consistency of Eqs. (fl"2"|) 
and (Ti3|) and can make the Arnoldi iteration to fail and 
produce complex energy solutions. 

The numerical errors in the matrix- vector products can 
be reduced by symmetrization. We thus calculate the 
RPA fields twice, first using the densities of a positive 
norm basis vector (X k ,y k ), and second using the den- 
sities of negative norm vector (y k * , X k * ) . The two re- 
sulting vectors are subtracted from each other to get the 
final stabilized RPA matrix- vector product, 
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to be used in Eq. ([12)) . Together, the Arnoldi iteration 
method and symmetrization of matrix-vector products 
stabilize our mean-field based iterative diagonalization. 



IV. TREATMENT OF SPURIOUS RPA MODES 

For the discussion of various spurious modes in the 
RPA method we refer the reader to, e.g., Ref. [l3j]. In the 
present study, we only consider spherical ground states 
neglecting pairing correlations, so the only spurious exci- 
tation is generated by the total linear momentum. There- 
fore, the only affected RPA mode is the isoscalar 1~ 
mode. In traditional RPA calculations that construct 
and diagonalize the full RPA matrix, the spurious 1~ 
mode is typically removed after the RPA diagonalization. 
Often a modified transition operator (|11|) is used, which 



has the property of (HF\ 



F, Pr, 



\HF) — 0, as long as 



the commutator is evaluated within a complete set of ba- 
sis states. In a finite model spaces of localized orbitals 
this relation is no more exactly valid, and the corrected 
operator does not remove spurious components exactly. 

To remove the spurious isoscalar 1~ mode from our 
physical RPA excitations we use the same method as 
in Ref. @, where the basis vectors are orthogonalized 
against the spurious translational mode P and its conju- 
gate "boost" operator R, which have the form: 
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The spurious RPA vectors {V,V*) and (11,11*) contain 
the p-h and h-p matrix elements of Eqs. ([20]) and (l2"Tj) . re- 
spectively. Our method differs from that of Ref. p] in the 
fact that we orthogonalize our basis during the Arnoldi 
iteration, which fits naturally with the iterative solution 
method and guarantees that the obtained approximate 
RPA excitations have exact zero overlaps with spurious 
modes. This is equivalent to diagonalizing the full RPA 
matrix in the subspace orthogonal to the spurious states. 
In our implementation, each generated new Arnoldi basis 
vector is orthogonalized as 
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where the overlaps A and fi are defined as 
(R, R*\X k ,Y k ) 



A 



(R,R*\P,P*) ' 
(P,P*\X k ,Y k ) 
~ (R,R*\P,P*) 



(23) 
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When more symmetries are broken, formulas equiva- 
lent to Eqs. (|22]) - ((2"4"]) can be used to remove spurious 
components coming from each broken symmetry of the 
mean field. 



CONVERGENCE OF STRENGTH 
FUNCTIONS 



The iterative Arnoldi method is meaningful for the cal- 
culation of strength functions only if the number of itera- 
tions needed for accurate results is significantly less than 
the full RPA dimension. To study how many Arnoldi it- 
erations we need for good accuracy, we calculated electro- 
magnetic isoscalar (IS) and isovector (IV) strength func- 
tions [IH for doubly magic nuclei. All calculations were 
performed by implementing the RPA iterative solutions 
within the computer program HOSPHE [hI, which solves 
the self-consistent equations in the spherical harmonic- 
oscillator (HO) basis. We studied both the convergence 
of smoothed strength functions as a function of number 
of Arnoldi iterations and as a function of the number of 
HO shells. 

We used the same definitions of the + , 1~, and 2 + 
transition operators as in Ref. [l7| and the Skyrme func- 
tional SLy4 of Ref. [20J . The function we used to smooth 
the strength functions was also the same as in [l6[ , with 
-Rbox = 20 fm. Because the HF ground state of 132 Sn 
is spherically symmetric, our approximate RPA phonons 
have good angular momentum. We tested the use of large 
basis sets up to 40 HO shells. The HF ground state ener- 
gies were well converged for all double magic nuclei when 
25 HO shells were used. Below, we present the results 
only for 132 Sn. 



Convergence as a function of the number of 
Arnoldi iterations 



1. The + strength functions 

Figure [Tj shows the + IS and IV smoothed strength 
functions for 132 Sn calculated with 100 Arnoldi iterations 
compared with the standard RPA results from Ref. [l7l |. 
Agreement between the strength functions is excellent. 
The four panels of Fig. [5] show the convergence of the 
smoothed strength functions of Fig. [TJ as the number of 
Arnoldi iterations increases. The panels show differences 
of the strength functions calculated at 20 iteration inter- 
vals. The convergence is quite satisfactory after 100-120 
iterations. 
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20 30 
E [MeV] 

FIG. 1: The + strength functions in 132 Sn calculated by 
using 25 HO shells and 100 Arnoldi iterations for the SLy4 
functional (solid lines) , compared with the standard RPA cal- 
culation of Ref. [13] obtained for the SkM* functional (dashed 
lines) . 
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FIG. 2: Convergence of the 132 Sn + strength functions of 
Fig. [T] Solid lines are for the IS and dashed lines are for 
the IV strength functions. Each panel shows the difference of 
two strength functions, one with n iterations and the other 
calculated with n — 20 iterations. 
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2. The 2 + strength functions 



3. The 1 strength functions 



Figures [3] and 0] show similar results as Figs. Q] and 
[21 but for the 2+ strength functions in 132 Sn. As for 
the + case, the IS and IV strength functions from the 
Arnoldi iteration agree very well with the strength func- 
tions of Ref. [l7[ • The convergence of strength functions 
is as fast as for + ; after 120 iterations the smoothed 
strength functions change only by about 5%. We thus 
have to make only 120 iterations to calculate reason- 
ably accurate 2 + strength functions for the RPA problem 
whose dimension is D = 1020. The large double spikes 
observed in Fig. [4] below lOMeV are due to the lowest 
RPA phonons, which by the smoothing procedure acquire 
~100-keV widths and move slightly down in excitation 
energy. 



Standard RPA 
Arnoldi n= 120 



: 200 




20 30 
E [MeV] 

FIG. 3: Similar to Fig. [Tjbut for the 2 + strength functions. 
All results were calculated for the SLy4 functional. 
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FIG. 4: Similar to Fig. [2]but for the 2 + strength functions. 



Figure [5] compares the 1~ strength functions of our it- 
erative method (solid line) with the strength functions 
from Ref. [TJj (dashed line). The solid line shows the 
IS strength calculated when the generated Arnoldi basis 
is orthogonalized against spurious mode during iteration 
and dotted line corresponds to similar iterative calcula- 
tion without orthogonalization. The low-lying state at 
0.72 MeV, which has a large overlap with the spurious IS 
1 _ mode disappears when the orthogonalization method 
of Eqs. (|22]) - (f2"3|) is used. Also for the 1" strength func- 
tion, 100-120 Arnoldi iterations were needed to produce 
reasonably accurate results, see Fig. [51 
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FIG. 5: Main panel: the 1~ strength functions in 132 Sn, 
calculated using 100 Arnoldi iterations and with spurious IS 
mode removed (solid line), and results of the standard RPA 
from Ref. [l7| ] (dashed line). Dotted line shows results of 140 
Arnoldi iterations without orthogonalization against the spu- 
rious IS mode. Inset: same as in the main panel, but for 
the IV strength functions. All results were calculated for the 
SLy4 functional. 



When no orthogonalization is made against the spu- 
rious mode, the obtained excitations contain small com- 
ponents of the spurious mode. This affects the physical 
part of the IS strength distribution, especially around 20- 
30 MeV. The standard RPA strength function of Ref. [H 
has not been corrected for the spuriosity but only the 
strength of the lowest-lying state that has a large over- 
lap with the spurious IS mode has been omitted. At 20- 
30 MeV, this strength agrees well with our uncorrected 
strength. 

The orthogonalization method improves the conver- 
gence of the strength function, because now the 8.3 MeV 
1~ excitation is lowest in energy and thus converges first. 
Without orthogonalization against the spurious mode, we 
need 140 iterations instead of 100 to get acceptably con- 
verged strength function. 
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FIG. 6: Similar to Fig. [4] but for the 1 strength functions. 
The IV strength-function differences were multiplied by the 
factor of 200 fm 4 . 
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FIG. 8: Similar as Fig. [4] but for the convergence of the 2 + 
strength functions as a function of the number of HO shells 
No. 



B. 



Convergence in function of the number of HO 
shells 



In Section E3 we showed our strength functions cal- 
culated with 25 HO shells. This was found to be satisfac- 
tory, and using more shells did not appreciably change 
the obtained strength functions. The only effect of using 
more oscillator shells was that we needed to use slightly 
more Arnoldi iterations to produce well converged re- 
sults. In the case of 40 shells, about twenty more itera- 
tions were needed, compared to calculations made with 
25 shells. In Figs. [3 [51 and[H we show the convergence of 
strength functions for the + , 2 + , and 1~ modes, respec- 
tively. Each panel shows the difference of two strength 
functions obtained in the intervals of ANq = 4 HO shells, 
between Nq of 22 and 38. 
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FIG. 7: Similar as Fig. [2] but for the convergence of the 
strength functions as a function of the number of HO shells 
No. 

These plots overstress the variations of strength func- 
tions in the sense that slight shifts of peaks create the 
oscillating patterns in the difference plots. To illustrate 
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FIG. 9: Similar as Fig. [6] but for the convergence of the 1 
strength functions as a function of the number of HO shells 
No- 



this point, in Fig. [10] we show the 1 strength functions 
calculated for N = 22, 26, 30, 34, and 38 HO shells. 
Poor convergence of the IS surface mode creates some 
uncertainty in the position and width of the high-energy 
bump. Larger bases should probably be used if converged 
results for this particular mode were required. 

As noted in Ref. [H , well before the maximum number 
of iterations (equal to the RPA dimension D) is reached, 
the iteratively generated RPA matrix in the Krylov space 
can become singular. In that case, the stabilized iteration 
method of Eqs. ((T2J)— (JTHJ) still protects us from obtain- 
ing complex RPA eigenvalues, but the condition number 
of the Krylov-space RPA matrix approaches infinity, be- 
cause one or more of its eigenvalues collapse nearly to 
zero. 

In the standard method, the RPA matrix is calculated 
by using the bare p-h basis states. In our method, we 
instead start from the pivot vector of Eq. ([2H)l and the 
Arnoldi iteration then produces the rest of our basis vec- 
tors composing the Krylov subspace. This subspace is 
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FIG. 10: Similar as Fig. [5] but for the 1~ strength functions 
calculated for the numbers of HO shells iVo = 22, 26, 30, 34, 
and 38 



TABLE I: Spherical RPA and QRPA dimensions D as func- 
tions of the number of HO shells No. 







RPA 






QRPA 




No 


0+ 


1" 


2+ 


0+ 


1" 


2+ 


10 


70 


195 


261 


390 


1040 


1510 


20 


205 


555 


766 


2880 


8180 


12720 


25 


273 


734 


1020 


5538 


15912 


25088 


30 


340 


915 


1271 


9470 


27420 


43630 



spanned by the eigenstates of the RPA matrix which has 
an overlap with the pivot vector. Thus, in general, the 
Arnoldi iterations can only be continued until this sub- 
space is exhausted in which case the condition number 
goes to infinity. However, with finite numerical preci- 
sion this maximum limit of Arnoldi iterations is further 
reduced. 

In a typical iteration, during the first few iterations 
the condition number of the Krylov-space RPA matrix 
fluctuates, then approaches a stable plateau, and finally 
suddenly goes toward infinity. When that happens, the 
iteration must be stopped and one must backtrack to 
the iteration where the condition number was still ac- 
ceptable. Therefore, the number of Arnoldi iterations 
can depend on the size of the HO basis, and the results 
presented in this section correspond to the numbers of 
iterations fixed according to this prescription. 



VI. SCALING OF ITERATIVE SOLUTION 
METHOD 



As can be seen in Table |H the RPA dimensions D 
for doubly magic spherical nuclei increase almost lin- 
early with the number of oscillator shells Nq. This is 
easy to understand, because in this case, only the num- 
ber of particle states increases and the number of hole 
states always stays constant. Therefore the time to solve 
the full RPA eigenproblem in this case scales approxi- 
mately as Nq. In the spherical QRPA, the dimensions 
scale roughly between N§ and Nq, and the full QRPA 
scales approximately as Nq or Nq. The physically inter- 
esting and computationally challenging calculations are 
for deformed nuclei with pairing, and we should therefore 
compare the A^o scaling of iterative and standard QRPA 
diagonalizations . 

In the case of all symmetries of the mean field being 
broken, the QRPA dimension D is: 

D = ^[(N + l)(N + 2)(N Q + 3)} 2 . (25) 

This dimension increases very steeply (Nq) as the num- 
ber of HO shells is increased. For No = 14, the QRPA 
dimension is D = 1849600, for example. The correspond- 
ing standard QRPA solution scales as A^ 8 and is thus 
untractable. Therefore, for the QRPA calculations in de- 
formed nuclei, we must truncate the single-particle space. 
The best method is to use the two-basis method [2l|, by 
which one solves the HFB equations in the basis gen- 
erated by the HF part of the HFB matrix, and truncate 
the basis using a cutoff on the obtained pseudo-HF single- 
particle energies. But even then, the QRPA calculations 
scale as the sixth power of the number of useful single- 
particle states, and are thus prohibitively difficult. 

500 



400 



-29.79+10.69N -0.877N^+0.034N^ 
-9.63+5.75H r 0.44N^-+O.017N^ 




FIG. 11: Times to calculate 100 Arnoldi iterations for the 
spherical QRPA method applied to 132 Sn as functions of No. 
Squares and circles show results for the 1~ and 2 + modes, 
respectively, and lines show cubic fits. 



We illustrate the benefits of the iterative solution of 
the RPA or QRPA equations over the traditional method 
by comparing how the numerical work increases in the 
iterative method as the HO basis is increased. 



In order to illustrate the scaling properties of the it- 
erative QRPA method, we calculated the corresponding 
matrix-vector products with our developmental QRPA 
code, where the pairing has been set to zero. The RPA 
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fields h only depend on the normal RPA density matrix 
and the calculation of pairing part of the QRPA matrix- 
vector product is very fast due to a small number of the 
pairing coupling constants and the simple density depen- 
dence of typical pairing EDF. Therefore the time to cal- 
culate the pairing part is negligible. The only significant 
increase of running time in spherical QRPA compared to 
RPA comes from the need to handle higher-dimensional 
basis vectors in the calculation of various overlaps and 
vector additions during the Arnoldi iteration. Therefore, 
as we keep all particle-particle and hole-hole RPA ampli- 
tudes in the calculation, but set them to zero, our tim- 
ing results accurately reflect the timing of the spherical 
QRPA calculation. 

In Fig.QTJ we show the scaling properties of the spher- 
ical QRPA calculation with iterative Arnoldi method. It 
is clear that the scaling of our iterative method is as Nq, 
that is, it is linear with respect to the QRPA dimen- 
sion D. Of course, the prefactor itself is linearly propor- 
tional to the number of Arnoldi iterations. However, as 
discussed in the previous section, the Arnoldi iteration 
method cannot in practice go full dimension before the 
generated Krylov-space matrices become singular. As 
long as we are satisfied with a few hundred iterations 
at most, the iterative method gives us a vast speed im- 
provement. In the full RPA or QRPA diagonalization, 
the calculation and storage of a very large dense RPA or 
QRPA matrices also takes a considerable additional time 
- a step that the iterative method avoids completely. 

In addition to the moment-method based iteration, 
which is ideal for strength functions, the iterative method 
can also be modified to be suitable for different kinds of 
other calculations. If we are interested in a number of 
very well converged lowest RPA eigenmodes, restarted 
Arnoldi methods 22j can be used. These methods use 
more iterations than basis states, i.e., after a maximum 
number d of basis vectors is generated, new approxima- 
tions for the wanted d' < d eigenmodes are calculated, 
and iteration is then continued to generate new improved 
basis states from d! + 1 to d again. The restarting can 
be made as many times as needed to produce wanted 
number of well converged lowest excitations. 

Methods such as Arnoldi or Lanczos produce conver- 
gence at the extreme ends of the excitation energy spec- 
trum. If eigenmodes away from the extremes are looked 
for, shift and invert methods [H, [24| can be used. These 
methods allow iterative methods to be used to find RPA 
eigenmodes anywhere inside the RPA excitation spec- 
trum. 



VII. SUMMARY AND CONCLUSIONS 

We have presented a method to calculate accurate RPA 
response functions by using the iterative Arnoldi diag- 
onalization related to the sum-rule conserving Lanczos 
method of Ref. [5J . We used strictly the same EDF for 
the ground state calculation and RPA excitations. We 
have showed how the Arnoldi method must be stabilized 
in order to apply it reliably to the RPA eigenvalue prob- 
lem. The resulting electromagnetic strength functions 
are in good agreement with the standard RPA results 
and are obtained with numerical effort smaller by orders 
of magnitude. 

Our method closely resembles the FAM of Nakatsukasa 
et al. @, 0] ! except that our iterative method is different 
and that we use the HO basis instead of the mesh in 
coordinate space. The FAM and our method both allow 
the existing EDF mean-field codes to be used for the 
calculation of the RPA or QRPA matrix- vector products. 
With minor modifications, mostly pertaining to the full 
implementation of the time-odd mean fields, these codes 
can easily be extended to RPA/QRPA. In particular, our 
future implementation of the deformed QRPA solution 
will be based on the code HFODD 

We also implemented the method to remove compo- 
nents of the spurious RPA modes from the calculated 
strength functions that keeps the physical excitations ex- 
actly orthogonal against the spurious excitations in any 
finite model space. 

The smaller numerical effort of the iterative Arnoldi 
method, and the fact that in this method one does not 
have to calculate and store the RPA or QRPA matrices, 
allows our method to be applied to the calculation of 
electromagnetic and beta decay strengths and strength 
functions for deformed heavy nuclei. Work to extend our 
formalism and codes to deformed superfluid nuclei is in 
progress. 
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